/* BoxMuller.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *
 */
#if !defined(_BoxMuller_)
#define _BoxMuller_

#include <src/algorithm/random/boost_sobol.h>
#include <src/math/norms.h>
#include <boost/lambda/lambda.hpp>

namespace boost {
namespace random {
/*! \brief Box Muller method.

	Fill a two dimensional vector according to the Box Muller method.
    Let \f$\{ x_i := \{x_i^k\}_{k=1,..,D} \in \RR^D\}_{i}\f$ be a sampling generator of the unit cube \f$[0,1]^D\f$ (default is a
    Sobol Generator).
    Rename and relabel this sequence, rejecting all sampling outside the unit ball :
    \f[
        \{ y_j := \{y_j^k = x_i^k\}_{k=1,..,D} \text{ if } \|x_i\|<1 \in \RR^D\}_{j}
    \f]
    This class is a generator described by the sequence \f$\{v_j := \{v^k\}_{k=1,..,D}\}_{j}\f$ filled with
    \f[
        v_j = y_{j} \sqrt{ \frac{-2s}{s}};s=\|y_j\|
    \f]
    Note : this method could be defined for higher dimensions.
*/
template <class T, class UniformGenerator = sobol_generator<32> >
 class BoxMuller : public generator_list<T,UniformGenerator>
 {
 public:
	 typedef generator_list<T,UniformGenerator> Base;
	BoxMuller() : Base(2) {};
	const T& operator()() // next sampling
	{
		OS_double s =2.;
		do {
			samples_.assign(-1.);
			samples_ += 2.*Base::operator()();
			s = jmmath::norm2(samples_);
		} while (s >1.0);
		if (s>0)
			s = std::sqrt(-2.*std::log(s)/s);
		else 
			s = 0.;
		samples_ *=s;
		return samples_;
	}
 };
}
}

#endif
